###########################################################################################
###### creating shape file versions of outputs ############################################
###########################################################################################

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(sp)
library(sf)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OutDirFin<-file.path(MyDir,"Outputs3_Final")
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_occs")
IndSppOutDir<-paste(OutDirFin, "/Final_Individual_Spp", sep="")
PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#set project name:
projectName<-"WHOSnakes"

#define selection criteria of interest for further analyses:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax[DatSum$Listing_comment!="not listed"]

#outcats<-c("Current","Median","Q90","Q10")
outcats<-c("Current")

worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))

#################################################################################
#make temporary temp dir for rasters in this job and set from default to this dir
dir.create(paste(MyDir,"/ShapeCombine",sep=""),showWarnings=FALSE)
rasterOptions(tmpdir=paste(MyDir,"/ShapeCombine",sep=""))
#################################################################################

for (cat in outcats){
  
  ###################################
  # general mapspp info #############
  ###################################
  
  for (sp in spp){ #Dendroaspis_polylepis still running
   
    print(paste(Sys.time(), "adding shapefiles for", sp, cat, sep=" "))
    
    n<-which(DatSum$gen_sp_subtax==sp) #gets row number for current species
    m<-which(spp==sp)
    mapsp<-DatSum$ModUnit[n]
    
    if(file.exists(paste(OccDir,"/",sp,"_allRecs.csv",sep=""))){
      sprecso<-read.table(paste(OccDir,"/",sp,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
      sprecs<-sprecso[sprecso$uncertainty<= (100000),]
      if(length(sprecs[,1])==0){
        sprecs<-sprecso
      }
      
      if(length(sprecs[,1])>0){     #only continue if there's any points
        
        #####################################################################################
        #Known areas of occurrence (buffers around occurrences cut to suitable habitat)
        
        KTOshp <- st_read(paste(IndSppOutDir,"/",sp,"_",cat,"_KTO.shp",sep=""))
        CDqsshp <- st_read(paste(IndSppOutDir,"/",sp,"_",cat,"_CDQs.shp",sep=""))
        SDMqsshp <- st_read(paste(IndSppOutDir,"/",sp,"_",cat,"_SDMQs.shp",sep=""))
        SDMcatshp <- st_read(paste(IndSppOutDir,"/",sp,"_",cat,"_SDMcat.shp",sep=""))
        
        
        
        #change modunit for each species... was wrong before:
        KTOshp$MODUNIT<-paste(gsub("_"," ",mapsp))
        CDqsshp$MODUNIT<-paste(gsub("_"," ",mapsp))
        SDMqsshp$MODUNIT<-paste(gsub("_"," ",mapsp))
        SDMcatshp$MODUNIT<-paste(gsub("_"," ",mapsp))
        
        
        
        
         Sys.sleep(0.1)
         print(paste(Sys.time(), "combining all polygons, adding", sp, sep=" "))
         
         #combine all shapefiles into one:
         if(m==1){
           
           st_crs(KTOshp)<-st_crs(worldmap)
           st_crs(CDqsshp)<-st_crs(worldmap)
           st_crs(SDMqsshp)<-st_crs(worldmap)
           st_crs(SDMcatshp)<-st_crs(worldmap)
           
           SDMfeature<-rbind(SDMcatshp,SDMqsshp,CDqsshp,KTOshp)
           
         } else {
           
           SDMfeature <- st_read(paste(IndSppOutDir,"/",cat,"_SDMfeatureAll_",m-1,".shp",sep=""))
           
           st_crs(KTOshp)<-st_crs(SDMfeature)
           st_crs(CDqsshp)<-st_crs(SDMfeature)
           st_crs(SDMqsshp)<-st_crs(SDMfeature)
           st_crs(SDMcatshp)<-st_crs(SDMfeature)
           
           SDMfeature<-rbind(SDMfeature,SDMcatshp,SDMqsshp,CDqsshp,KTOshp)
           
         }
         
         print(paste(Sys.time(), "saving", sep=" "))
         
         st_write(SDMfeature, paste(IndSppOutDir,"/",cat,"_SDMfeatureAll_",m,".shp",sep=""),delete_layer=TRUE)
         file.remove(paste(IndSppOutDir,"/",cat,"_SDMfeatureAll_",m-2,".shp",sep=""))
         file.remove(paste(IndSppOutDir,"/",cat,"_SDMfeatureAll_",m-2,".dbf",sep=""))
         file.remove(paste(IndSppOutDir,"/",cat,"_SDMfeatureAll_",m-2,".prj",sep=""))
         file.remove(paste(IndSppOutDir,"/",cat,"_SDMfeatureAll_",m-2,".shx",sep=""))
      }
    }
  }
}

st_write(SDMfeature, paste(IndSppOutDir,"/",cat,"_SDMfeatureAll.shp",sep=""),delete_layer=TRUE)

######################################################
#delete the temp dir for this scriptL
unlink(paste(MyDir,"/ShapeCombine",sep=""), recursive = TRUE) 
######################################################

# #########################################
# # combining polygons into one shapefile #
# 
